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We suggest a new approach for transport through finite systems based on the Liouville equation. 
By working in a basis of many-particle states for the finite system, Coulomb interactions are taken 
fully into account and correlated transitions by up to two different contact states are included. This 
latter extends standard rate equation models by including level-broadening effects. The main result 
of the paper is a general expression for the elements of the density matrix of the finite size system, 
which can be applied whenever the eigenstates and the couplings to the leads are known. The 
approach works for arbitrary bias and for temperatures above the Kondo temperature. We apply 
the approach to standard models and good agreement with other methods in their respective regime 
of validity is found. 



I. INTRODUCTION 

Transport through nanosystems such as quantum dots 
and molecules has received enormous interest within the 
last decade 1,2,3 . Typically this problem is treated within 
one of two different approximations: (i) Rate equations^ 
for electrons entering and leaving the system, which can 
also take into account complex many-particle states in 
the central region^. Here broadening effects of the lev- 
els are entirely neglected. It can be shown that these 
rate equations become exact in the limit of high bias 7 , 
(ii) The transmission formalism, which is usually evalu- 
ated by Green function techniques^ (alternatively, scat- 
tering states can be calculated direct!}^) , allows for a 
consistent treatment of level broadening due to the cou- 
pling to the contacts. In principle, many-particle ef- 
fects can be incorporated into this formalism, but the 
determination of the appropriate self-energies is a diffi- 
cult task, where no general scheme has been found by 
now. Thus, many-particle effects are usually considered 
on a mean-field basis including exchange and correlation 
potentialsii*i2iiii4, which are of particular importance 
for the transport through molecules. Mean- field calcu- 
lations are well justified for extended systems, such as 
double-barrier tunneling diodesi^^, which exhibit many 
degrees of freedom (e.g., in the plane perpendicular to the 
transport). However, the bistability frequently obtained 
for such structures is questionable for systems with very 
few degrees of freedom as studied here. See, e.g., the 
discussion in Sec. III.B.4 of Ref. HH 

In our paper we want to bridge the gap between these 
approaches by considering the Liouville equation for the 
dynamics of the central region coupled to the contacts. 
The approach works within a basis of arbitrary many- 
particle states, thus fully taking into account the inter- 
actions within the central region. While the first or- 
der in the coupling reproduces previous work using rate 
equations^, the second order consistently takes into ac- 
count broadening effects. This is analogous to the con- 
sistent treatment of broadening for tunneling resonances 
in density-matrix theory^ 

The paper is organized as follows: We first present the 
formalism in section^] Then we demonstrate its applica- 



tion to the simple problem of tunneling through a single 
level, section lLTTI We show explicitly that the exact Green 
function result is recovered for all biases and tempera- 
tures. In sections llVl we give results for the double-dot 
system with Coulomb interaction where both standard 
approaches fail. Finally we consider the spin-degenerate 
single dot in section to investigate Coulomb blockade 
as well as the limit of low temperatures. 

II. INTRODUCING THE FORMALISM 

The total Hamiltonian for the system consisting of 
leads and the dot can be written as 

H = Hvy + -^Lcads + Ht- (1) 

The first term describes the dot. Our key issue is the 
assumption that the dot can be diagonalized in ab- 
sence of coupling, and the (many-particle) eigenstates 
and eigenenergies for Hp are denoted \a) and E a . Thus 
we have 

H D = Y J E a \a){a\. (2) 

a 

The leads are described by free-particle states 

#Lcads = ^ E kcTlC kul C k(Tl (3) 

where a =f, | describes the spin, k labels the spatial wave 
functions of the contact states and £ denotes the lead. In 
the following we assume two leads, i.e. £ = L,R, but 
generalization to more leads is straightforward. Finally, 
the last part in the Hamiltonian expresses the tunneling 
between the states in the leads and the dot 

H t= E [T b a(ka£)\b)(a\c kae + cl^b\T b * a (ka£) . 

kcr£,ab 

The matrix element Tb a (ka£) is the scattering amplitude 
for an electron in the state ka£ tunneling from the lead 
onto the dot, thereby changing the dot state from state 
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| a) to a state \b). Their evaluation is sketched in Add. 1X1 
Note that this amplitude vanishes unless the number of 
electrons in state \b), Nb, equals N a + 1. We will gener- 
ally denote states such that the particle number increases 
with the position in the alphabet of the denoting letter. 

Before proceeding it is important to introduce a consis- 
tent notation in order to keep track of the many-particle 
states in the leads. A general state vector for the entire 
system is written as \ag) = \a) ® \g), with \g) — \{N k e a }) 
denoting the state of both leads where Nkat € {0, 1}. 
Throughout the derivation of the general equations we 
use the following notation to ensure the anti-commutator 
rules of the operators 

• \g - ka£) = c kal \g) and \g + kal) = c{ at \g). 

I.e. \g — kai) denotes the same set of indices as 
the state \g), but with N ka -i reduced by one. Fur- 
thermore it contains a minus sign depending on the 
number of occupied states to the left of the position 
kai. 

• \gkai) = c\ al c kal \g) and \gkal ) = c kal c\^\g). 
I.e., \gkai) = S Nk<reA \g). 

• The order of indices is opposite to the order of the 

operators. E.g. \g - k'a't' + kal) = c%<rt c k' a' = 

-c kla ,A*l\9) = ~\9 + kai - k'a't) for kai + 
k'a'l 1 ', which is tacitly assumed, unless stated oth- 
erwise. 

To simplify the notation, at is only attached to k the 
first time the index k appears in the equation, and in the 
following it is implicitly assumed to be connected with k. 
We also use the convention that ^2 kcr m means summing 
over k and a with a fixed i, which is being connected to 
k in this sum. 

The matrix elements of the density operator p are de- 
noted Pag-bgi = ( a g\p\bg') an( i the time evolution of the 
matrix elements are governed by the von Neumann equa- 
tion 

m ^ t Pag,bg' = ( a 9\ H P ~ P H \ h 9') ( 5 ) 

The particle current from the left lead into the struc- 
ture, Jl, equals the rate of change in the occupation of 
the left lead. We find that 

J t = ~i E ( c K) = £ pw 

kcr(L) ka(L) 

2 f 1 (6) 

ka{L),cb Kg ) 



where we have used the definition of the density operator 
to calculate the average value of the number operator in 
the left lead. 

The goal is to determine these elements of the density 
matrix, which describe the correlations between the leads 
and the dot. They are determined using the equation-of- 
motion technique, and from Eq. J5J we obtain 

ih-T^Pcg-kviibg = ( E c ~ E b - E k )p cg _ k . bg 

+ E ^ cb ' (k)Pb>gk:bg — E Pcg-k;c' g-k^c'bik) 
V c> 

+ E [j2 T cb'( k ')Pb'g-k+k';bg+J2 T dc( k ')Pdg-k^k';bg 
k'a'i' b' d 

~ E Pcg-k;c' g-k'Tc'b(k') — Pcg-k;ag+k'Tb*a(k') . 
c' a 

(7) 

While p cg _ klJ £.bg describes the transition of an electron 
with quantum number k and spin a from lead I to 
the central region, terms like Pb> g -k+k'-bg describe the 
correlated transition of two electrons with k and k' . 
Pb'g-k+k'-bg satisfies a similar equation of motion con- 
taining also correlated transition of three electrons on 
the right hand side. In order to break the hierarchy we 
apply three approximations: 

(i) We only consider coherent processes involving tran- 
sitions of at most two different fc-states. (ii) The time 
dependence of terms generating two-electron transition 
processes is neglected, which corresponds to the Markov 
limit. 20 (iii) We assume that the level occupations fkal in 
the leads are unaffected by the kinetics of the dot, so it is 
possible to factorize the density in the leads and on the 
dot. This is realistic for "large" leads which are strongly 
coupled to reservoirs, i.e. good contacts. 

Defining 

W b'b = E^'ffibff' ( 8 ) 
9 

^ba(kai) = J2Pb9-k;ag (9) 
9 

we find the following set of coupled differential equations 
(see Appendix for a detailed derivation): 
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ih—(f> cb (ka£) =(E C - E b - E k )4> cb (k) + ^^T cb >(k)f k w b > b - ^(1 - fk)w cc >T c , b (k) 



E 



a,b' ,k' a' £' 



E 



a.b' ,k' a' £' 



E 

a,b' ,k' a 1 'I' 

E 

h'c'k'a'l' 



E 

b',c',k'<r'i 

E 

c' ,d,k' a' 

E 

c',d,k'a'i' 

E 

c' ,d,k' a' £' 



\T cb .(k')j k ,<\> bla [k) - (1 - j k )U\k')T b , a [k)\ T* a (k>) 
E k + E k , - [E c - E a ) + iO+ 

[(1 - .fk^c b '{k)T b , a {k') - T cb ,(k)f k b/a (k<)} T b * a (k') 
E k + E k , - (E c - E a ) + iO+ 

T cb> {k') lf k ^ b ,g(k)T b * a (k>) - T b>a (k)f k r ba (k')} 
E k - E k > - (E b , -E b )+ iO+ 

T cb >(k!) [T c *, fc ,(fc')(l - fj^m - (1 - hW c , b ,(k')T c , b {k)] 
E k - E k > - (E b , -E h ) + iO+ 

[fk4cv(k)T* clb ,(k>) - T cb ,{k)f k r c , b ,(k')] T db {k') 
E k -E k , - (E c -E c ,) + iO+ 

[T: c (k')(l - f k >)(pdc>{k) - (1 - f k )^ dc {k')T dc ,{k)]T c , b {k') 
E k -E k ,-{E c -E c ,)+m+ 

T2c(k') [(1 - fk-)<j>dc>{k)T c , b {k') - T dc ,(k)f k <p c , b {k')] 
E k + E k , -{E d -E b ) + iO+ 

TUk') [T dd {k')f k ,4> db {k) - (1 - f k )4> dcl {k')T db {k)] 
E k + E k , - [E d - E b ) + iO+ 



(10) 



iTi—ww = (E b - E b >)w bb > 



■ E (T ba {k)4>l, a {k) 

a,kcr£ 



4> ba (k)T b *, a (k)) + J2 (T: b (k)4> cbl (k) - 4>l b {k)T cV {k)) (11) 



c,ka£ 



These equations are the main result of this paper. They 
satisfy current conservation, as shown in App. El The 
numerical implementation of this approach is straightfor- 
ward and we will give examples in the following sections. 

If we entirely neglect the correlated two-particle tran- 
sitions, only the first line of Eq. I|10[l remains. Applying 
the Markov limit we obtain a set of equations analogously 
to Eqs. (2a, b) of Ref.0 This shows that our approxima- 
tion scheme goes substantially beyond the rate equation 
scheme of Gurvitai^, which only holds in the high-bias 
limit. 



Inserting this in Eq. I|10|) with c = 1 and b = gives 



ih— 4> W (M) = [E x -E k + Z(E k )} ^o(H) 



- E k -E' k +W 



Tffl £ ^ZrVJ + Ti{k){h - tan), 



k't 



where the self-energy 



\Ti{k)\' 



kf 



E-E k + i0+ 



(12) 



(13) 



III. SINGLE LEVEL WITHOUT SPIN 

In order to demonstrate the formalism described in the 
previous section we consider a single level without spin. 
We show that this case can be solved analytically in the 
stationary state and that the exact nonequilibrium Green 
function result is recovered. 

The possible dot states are the empty state with 
energy E — and the occupied state 1 with energy E 1 . 
The coupling matrix elements between the leads and the 
dot are Tio(k£) — Ti(k), and the others equal zero. 



has been introduced, and we have used the normalization 
of the probability wqo + wn = 1. 

After multiplying Eq. (jf^ with T^(k)S(E - E k ) and 
summing over all fc-states (in a fixed lead £) we obtain 



ih-B( (E) = [E X -E+ £(£)] B( (E) 

Te(E) f dEl Bft(E') + B$(E') (w) 



+ 



2tt 
Tt(E) 
2tt 



E-E> + i0 H 
[ft(E) - w n ] 
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for the new variable 

Bi Q (E) = J25(E-E k )TZ(k)fao(M) (15) 

where T e (E) = 2nJ2 k 8{E - E k )\T e (k)\ 2 . 
Eq. 111(1 becomes 

=-~ / AE^{B^{E) + B^{E)} . (16) 

Finally, the current formula Eq. @ yields 

Jl = ~lj dE ^ {B^(E)} . (17) 

Throughout this paper, we apply Fermi functions 
hat = l/[exp((E k ~ li t )/k B T) + 1] = f t (E k ) for the 
lead occupations with chemical potentials \xg and tem- 
perature T. Except for this section, the bias V is applied 
symmetrically around zero, i.e. (i^ — V/2, fip/ = — V/2. 
The contact functions Ti(E) are assumed to be zero for 
\E\ > W , while they take the constant values 1^, inde- 
pendent of spin, for \E\ < 0.95VF. For 0.95W < \E\ < W 
we interpolate with an elliptic behavior in order to avoid 
discontinuities. 

The time-dependent net-current Jr (t) flowing from the 
right lead into the single level has been calculated from 
Eqs. I|14lltjll7|) in the following situation: For times t < 
the chemical potentials of both leads and the single level 
are aligned, i.e. fi° L = u R = E\ = 0. At t = the chemi- 
cal potential of the left lead is raised instantaneously to 
giving a step- like modulation of the bias. The re- 
sult is shown for different values of fiL in Fig. ^ Also 
shown is the result of an exact time-dependent Green 
function calculation^ It is not surprising that our re- 
sults do not show the exact time-dependence because the 
Markov limit has been invoked in the derivation of the 
generalized equation system in Eqs. (110111(1 . 

In the long-time limit, we reach a stationary state with 
the current 

T = T =- f dE T L{E)T R {E)[f L {E) - f R (E)} 
L R hj 2tt l^-^ -£(£)| 2 

(18) 

which is derived analytically in App.[Dj Eq. ((18(1 is in full 
agreement with the exact nonequilibrium Green function 
result. 8 
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FIG. 1: (Color online) The time-dependent current calcu- 
lated with our method (full line) and with the time-dependent 
Green function method from Ref. (dashed line) as response 
to a step-like modulation of the bias with step height fiL- The 
coupling is Tl = Tr = T/2, the temperature ksT = 0.05r, 
and the half- width of the band is W = 30I\ 



the electron by a high magnetic field.) Denoting the 
left/right dot by a//3, the Hamiltonian reads 

H =E a S a d a + Eprfpdp + Udld a d^d 

+ Y {}kL d i c kL + hR d l c kR + ft-C-) 
k 

with f2 being the interdot tunneling coupling and U the 
Coulomb energy for occupying both dots. The first four 
terms describe the isolated double quantum dot Hjj . Di- 
agonalizing this part of the Hamiltonian gives the follow- 
ing states 





= 10), 


E 


= 0, 






11) 


= (ai<& + /3idJ)|0>, 


Ei 


1 

~ 2 


P 


- VA 2 + 4^ 2 


|2> 


= (a 2 di+(3 2 dl)\0), 


E 2 


1 

~ 2 


P 


+■ VA 2 + 4ft 2 


\d) 




E d 


= E a 


+ 


Ep + U, 



IV. DOUBLE QUANTUM DOT 

The double quantum dot structure, where the dots are 
coupled in series, is a standard example to study tunnel- 
ing through a multiple-level system. In case of Coulomb 
interaction and finite bias the validity of both the rate 
equation method and the Green function formalism is 
limited. 

To simplify the analysis we treat the spinless case. (A 
possible realization is to favor one spin polarization of 



with A 

"1/2 = 



E n 



-Ep, P = 
and Pi/2 



E a + Eg, C± 



A±VA 2 +4Q^ 
20. ■ 



From Appendix 1X1 we find 



Hr — 



[Tl (k)\l)(0\c u +T« (k)\2)(0\c H 



T l dl (k)\d){l\c M +T l d2 (k)\d){2\c M 



(20) 
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with (skipping the k-dependence of the matrix elements) 



-L 1 1 



^(i|4lo> 



t L a t , 



rpL 



t* L (2\dl\0)=tla* 2 



= *I<d|4|l) = tI/3 1 , 
= *«<d|4|l> 



r& = *£<d|4|2> = tlfc, 

T& = t* R {d\d\\2) = -t R a 2 , 



where the signs of the coupling matrix elements are due 
to the order of the operators in the double-occupied state. 

Applying the method in the same way as for the sin- 
gle level system gives eight different functions of the type 
B~AE) and five different occupations Ubf- These equa- 
tions have been solved and the stationary current has 
been recorded. 

By comparing with exact Green function results, it 
has been verified numerically that in the non-interacting 
case (U = 0) the exact transmission is obtained for var- 
ious values of level splitting and interdot coupling (not 
shown). Furthermore, for both U — and non-zero U 
we have calculated the stationary current in a situation 
where the levels are de-aligned with E a — —Ep = 0.5T, 
and f2 = I\ The results for different values of U are 
shown in Fig. [5] together with the Green function result 
for U = 0. Obviously, the latter is fully recovered in 
the non- interacting limit. The straight dashed line in 
the figure is the quantum rate equation result obtained 
by Stoof and NazaroA^, which is valid in the high-bias 
limit (V — > oo) for U — > oo. The same result is found 
in Ref. 7 using another rate equation method. The small 
discrepancy between the results could be due to the finite 
bandwidth used in our calculation. For intermediate val- 
ues of U the results looks reasonable and exhibit a smooth 
interpolation between the limiting cases. The kink on the 
curve for finite U is due to the single occupied state. 



V. SPIN-DEGENERATE LEVEL 

Now we consider a spin-degenerate single level with 
energy E\ and Coulomb interaction U. We use the pa- 
rameters U = 1.9 meV, and T = T L + T R = 0.295 meV, 
as experimentally determined for the structure studied 
m Ref.lH The conductance 



G = e 2 



dJ 



d(^L - mr) 



(21) 



is expected to reach G — \ 8T L T R /(T L + T R ) 2 in the 
zero-bias limit /i^ — > fi R for temperatures far below the 
Kondo temperature T~a 24 ' 25 . As G max w 0.5e 2 /h in the 
experiment we use Tl = 0.275 meV and T R = 0.02 meV. 
Furthermore the band width W — 5 meV is applied. 
In Fig. |3 we show the zero-bias conductance as a func- 
tion of the dot level, which is modified by a gate bias 
in the experiment. We find the standard Coulomb os- 
cillations, where the conductance exhibits peaks when- 
ever the single-particle excitation energies are close to 



0.2 



0.15 



0.1 



0.05 - 



1 












- U/r = (Green tunc.) 






— u/r = o 






— u/r = 2 




II 


u/r = 5 






u/r= 15 






Stoof and Nazarov 






i 





10 



v.. /r 

bias 



15 



FIG. 2: (Color online) Stationary current through the dou- 
ble quantum dot structure for different values of the interdot 
Coulomb repulsion U. The triangles are from a nonequilib- 
rium Green function calculation, and the dotted line is the 
result by Stoof and Nazarov— valid in high-bias limit for 
U — * oo. The levels of the dot are placed symmetrically 
around the zero- bias with E a — Ep = T. We use the interdot 
tunneling coupling Q, — V, Tl — T_r = T/2, the temperature 
k B T = O.ir, and the half-width of the band W = 20I\ 




FIG. 3: (Color online) Zero-bias conductance as a function of 
level position for different temperatures. All parameters are 
according to the experimental data shown in Fig. 2 of Ref . I23I 



the Fermi edge of the contacts, \i — (depicted by ver- 
tical dashed lines at E\ = and E\ = —U). The peak 
positions and widths are in good agreement with the data 
given in Fig. 2 of Ref. l23l The peak heights for the peak 
around E\ « agree reasonably with the experiment, if 
one takes into account that for elevated temperatures the 
presence of different levels raise the conductance which is 
not included in our single-level model. (The experimen- 
tal level spacing corresponds to 5 K.) The experimental 
peak heights for the peak at E\ as — U are lower, while 
they are exactly identical with the corresponding peaks 
E\ w due to electron-hole symmetry in our calcula- 
tion. Possible sources for this deviation result from an 
energy-dependence of the Te(E) in the experiment or the 
admixture of different levels. 
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Ml-^R (meV) 

FIG. 4: (Color online) Differential conductance for finite bias. 
Parameters as in Fig. [3] 



Further lowering the temperature, the zero-bias con- 
ductance should increase in the region < E\ < U, due 
to the Kondo eSed&&. Albeit we observe an increase 
in parts of this region, the (probable unphysical) dip in 
our curve for T = 0.1 K at E\ = —{7/2 persists even at 
lower temperatures. Furthermore the conductance can 
exceed Go at the peaks. This indicates that our ap- 
proach fails in the Kondo limit, where strong correlations 
between lead and dot state require elaborated renormal- 
ization group2£i2L2& or slave bosonSLSI techniques. 

In Fig. 0| we show the finite bias conductance at 0.8 K, 
where both the conductance peaks for fi^ ~ (J-r discussed 
above as well as the excitations can be detected. We ob- 
serve a strong asymmetry due to Tl T^j. This can 
be understood from Fig. [S] where the current is plotted 
versus bias at the single-electron excitation peak E\ = 0. 
For negative bias the electrons rapidly leave the dot via 
the thin left barrier and the dot is essentially empty. 
Thus both spin directions can tunnel through the thick 
right barrier, which is limiting the current. In contrast, 
for positive bias the dot is occupied with a single electron 
(as long as hl — V/2 < E\ + U) with a given spin and 
only this spin direction may tunnel through the thick bar- 
rier, reducing the current approximately by a factor of 2. 
We have shown the respective results for the rate equa- 
tion modeli for comparison. The short-dashed horizontal 
lines refer to a bias which allows only single occupation 
of the dot, while the long-dashed line considers the case 
where both the single- and the two-particle state are lo- 
cated between both Fermi levels. The currents from the 
rate equation model slightly exceed our results, as the 
peaks are not completely within the bias window due to 
broadening. 



VI. DISCUSSION AND SUMMARY 

We have presented an approach for transport through 
finite systems based on the Liouville equation. This 
approach recovers the results from the Green-function 
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FIG. 5: (Color online) Current versus bias for E\ = and 
T — 0.8 K. Parameters as in Fig. [3] The dashed horizontal 
lines correspond to the rate equation model by Gurvitz and 
Pragej-, where we added the respective formulae. 



method in the noninteracting limit for the models stud- 
ied. In the high-bias limit the results are consistent with 
the many-particle rate equations. Thus it bridges the 
gap between these approaches and allows for a consistent 
treatment of Coulomb interaction and broadening effects 
for arbitrary bias. E.g. Coulomb blockade peaks are 
correctly reproduced. The model fails below the Kondo 
temperature where strong correlations between the finite 
system and the contacts dominate the behavior. 

Correlations between tunneling events have been previ- 
ously studied by the method of a diagrammatic real-time 
technique^. While this work was completed we also be- 
came aware of a cumulant expansion of the tunneling 
Hamiltonian 32 . It would be interesting to study the rela- 
tion between these approaches and our method. A cen- 
tral question is here, whether the exact Green function 
result, such as Eq. JTSJ, can be obtained for the nonin- 
teracting case. 

The numerical implementation of our approach is 
straightforward and explicit results were presented for 
standard model systems made by up to two single- 
particle states. For larger systems the number of many- 
particle states 6, c increases dramatically, and so does the 
number of <^b c -functions. Thus sophisticated routines are 
needed for the implementation and evaluation of real sys- 
tems. 
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APPENDIX A: DETERMINATION OF MATRIX 
ELEMENTS T ab (k) 

Conventually one starts with a single-particle basis in 
the central region with wave functions ^(r), spin func- 
tions Xrr and associated creation operators S n(J . Then an 
arbitrary many-particle state \a) can be written as 

i«> = E fl j44---4vJ > 
j 

where ji — rt{Oi determines the i-th single-particle state 
in the iV a -particle Slater determinant determined by the 
index set j = (ji , ji , ■ ■ ■ JN a ) ■ In order to avoid double 
counting, we restrict to the ordering ni < n% < . . . < 
riN a , where spin-up precedes spin-down for equal n. The 
expansion coefficients aj can be obtained by exact diag- 
onalization of the dot Hamiltonian. 

In the single-particle basis the tunneling Hamiltonian 
reads 

H T = E \Xii{k(rt)c\ al d n(T + tniikcrfydl^Ckvi] (Al) 

kcr&,n 

Inserting the unit operators ^ a l a )( a L I^X^I we nn d 



kcri.a.b 



=T* fa ( kcri 



i)J2tne(ka£)(a\dt a \b)(b\ 



Ckal 



(A2) 



=T ab (ka£) 



to be used in Eq. Q. 



APPENDIX B: DERIVATION OF EQS. (1101111) 

Using the approximation (i) we find for one of the two- 
electron transition terms in Eq. Q) 



~d? 



^ — Pb'g-kai+k'a'l'-bg = i E b' + E k' ~E b ~E k )Pb'g-k+k' -b 



~ E T b'a (k)Pagk+k> -bg + E T o'V ( fc ') P d g-kV-bg 
a c' 

"Y! ft'a- Wlc'a-t^'ii W~y! Pb' g-k+k' ;ag+k ,r Eba{k') 



Now we take the Markov limit (ii) following the stan- 
dard treatment of density matrix theory for ultrafast 
dynamic^. This implies adding — iO+p t)/g _ jfc(r£+fc , cr ,^. { , g 
on the right hand side in order to guaranty the decay of 
initial conditions at t — — oo, and neglecting the time de- 
pendence of the inhomogeneityS (second and third line) . 



Then this linear differential equation can be solved di- 
rectly, yielding 



Pb'g-kcrl+k'o't'-.bg 



1 



E k - E k , - {E v - E b ) + 10+ 



X [ -J2 Tb ' a WP*9k+k<;bg +J2 T *'b>( k ')P c >g-kW:bg 
a c' 

Pb' g-k+k' :d g-kEc'b{k)~^2 Pb' g-k+k' :ag+k'Tba{k') 



In the same way the other two-electron transition terms 
in Eq. (J7J are determined by 



Pdg-kcrl-k'cr't'\bg 



E k + E k , - {E d - Et) + i0+ 
: E [ — Edc'(k)p c i gk _ k /. bg + Tdc'{k')p c ig_ kk i.^ g 

c' 

Pdg-k-k' ;c' q-kEc'b{k) — Pdg-k-k' ;c' g-k'Tc'b(k') 



Pcg-ka£;c'g-k'a't ^ _ Ey _ ^ _ £^ + iQ+ 

X [ E ^c&' (tyPVaklc'Q-k' + E Tdc{k')Pdq-k-k';c'q-k' 
V d 

~ E Pcg-k-dg-k'-kEdc'(k) — E Pcq-k-b' gk'^c'b'jk') 



Pcg—ka£\ag+k'a't' 



E k + E k , - {E c - E a ) + i0+ 

E \Ec b '{k)p b , gk . ag+k , + T cb '{k')p b , g _ k+k ,. ag+k , 
V 

- Pcg-k;b'g+k'-k T b'a{k) - P ca _ k , b , ak T T b' aik') 



In order to obtain Eq. (|10|) we sum over g in Eq. (JJJ 
after inserting the above approximations for the two- 
electron transition terms. Using the definitions (|8I9|) and 
the decoupling assumption (iii) we obtain 



^2pb'gk;bg = E^' 1 



Pb'g-.bg 



fk E Pb'g-bg = !kW b > b 



Similarly 



J2Pb'g k ;bg~ (I" fk)w b 'b, 

g 

^2jPbg-k'k;ag ~ fk<f>ba(k'), 

g 

E Pbg-k'k;ag ~ ( l ~ fk)<Pba{k') 
9 

Furthermore note that 

^^Pbg-.ag+k = E Pbg-k;ag = <t>ba(k) 
9 9 

E P b ' '9+k;bg+k = E ^'3 fc ; b 9 ~ IkWb'b 
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as well as similar relations hold, where \g) is identical 
with \g) except for exchanging 1 and in the occu- 
pation of state k (including the appropriate change of 
sign). Particular care has to be taken in order to insure 
the anti-commutation rules. E.g., Y. g Pb< g -k+k>:ag+k> ~ 

~fk'4>b'a{k)- 

In the same way 

'^-JlPbg-b' g = ( E b - E b')Pbg-b'g 

+ Tb a (k)p ag+k . b , g + 2^ ^cb(^)Pcg-k;b'g 

a,ka£ c,ka£ 

Pbg;cg-kTcb>(k) ~ Pbg-ag+k^b* a(^) '• 

I ajzcji 

gives Eq. after summing over g. 

APPENDIX C: CONSERVATION OF CURRENT 

We will in the following show that the formalism obeys 
current conservation, i.e. 

(CI) 



c,ka£ 



-£{N D ) = J L + J R , 



with Nj) = ^2 b Nt,\b)(b\ being the number operator of the 
dot. From the definition of the density operator we get 



pN D ] =Y, N »at 



w bb . 



(C2) 



The time derivative of Wbb is obtained from Eq. 1 (11(1 



df 



Wbb 



akat 



^3{T c * fc (fc)0 cb (fc)} 



(C3) 



cka£ 



Inserting this in Eq. (|C2|) and renaming the summation 
indices in the second term leads to 

sW=lE^- N a p{T b * a (k)<t> ba {k)} . (C4) 

bakal 



Now the T ba (fc)-matrix elements are vanishing for Nb ^ 
N a + 1, and the right-hand side of Eq. l|C4Jl becomes 
Jl + Jr using the definition of the currents Eq. Thus, 
current conservation 1C1|) holds. 

APPENDIX D: DERIVATION OF EQ. fTgf) 

Defining B 1Q = Sf + B? and V = T L + T R we find 
from Eq. ifHfl 



ih—B 10 (E) = [E 1 -E + $t{X(E)}]B 10 (E) 
+ T(E)Z {B 10 (E)} - ^ / dE' V I B * w{E ' 



2tt 



T L (E)f L (E)+T R (E)f R (E) 
2n 



E-E' 

r(E) 



2tt 



(Dl) 



where 3 {£(£;)} = -Y{E)/2 has been used. Eq. (JdTJ is 
a linear inhomogeneous differential equation which has a 
particular stationary real solution i3fp at (i?) determined 
by 



m r , f Bf at (E') \ 

2n J \ E-E' J 



= [E 1 — E + {£(£)}] Bf^iE) 
T L (E)f L (E)+T R (E)f R (E) T{E) 



2n 



2tt 



(D2) 



Numerically, we find that this solution is indeed reached 
from different initial conditions in the long-time limit. 
Inserting the integral over Bf ( ^ t (E) from Eq. 1D2I) into 
Eq. Ijl 4(1 gives the stationary solution 



[E.-E + ^E^B^E) 
= r L (E) 



E X -E + V {E) M 
T(E) 



Bf at {E) 



T L {E)T R (E)[f R {E)- f L (E)] 
2t:T{E) 



(D3) 



As Bf 1 ^ (E) is real it does not contribute to the imaginary 
part of B^ (E) in Eq. IT7|l providing the final result {T%)l . 
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